Journal of Theoretical Biology

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Moment-closure methods provide a widely used approach for approximating low-order statistics of stochastic biochemical reaction networks, but the resulting closed equations may exhibit physically meaningless dynamics. Here, we systematically investigate the stability of second-order moment-closure approximations for a bursty gene expression model with nonlinear protein degradation. We consider the Gaussian, log-normal, and negative-binomial closures. By analyzing the resulting dynamical systems, we show that the Gaussian closure fails the proposed stability criterion for all parameter values, whereas both the log-normal and negative-binomial closures possess a unique physically accessible fixed point and are stable if and only if kB(1 + 2B) [≥] 1, where B is the mean burst size and k is the dimensionless protein synthesis rate. Moreover, this condition is automatically satisfied when the steady-state mean protein copy number is at least one, suggesting that it is mild in biologically relevant regimes. Numerical comparisons over a broad parameter range further show that the negative-binomial closure is the most accurate, while the Gaussian closure is the least accurate. We further develop higher-order negative-binomial moment-closure schemes and demonstrate that increasing the closure order substantially improves predictive accuracy. These results highlight the importance of considering dynamical stability, in addition to numerical accuracy, when assessing moment-closure approximations.

Fractional differential equations are increasingly used to represent temporal memory in tumour immune and pharmacokinetic models, but changing the derivative order without a consistent time scale can make comparisons difficult to interpret. We study a dimensionless, commensurate Caputo system for tumour, immune, normal cell, and drug states. We establish positivity, boundedness, the tumour invasion threshold, and local stability conditions. We then use a fractional Adams Bashforth Moulton predictor corrector method with an explicit mesh convergence study. Two conclusions correct potentially misleading interpretations of the model. First, the equilibrium locations are independent of the fractional order: decreasing the order changes transient time scales but does not by itself alter the model's asymptotic states. Second, after a unit bolus the drug equation has the exact solution C({tau})=E(-K{tau}). Its maximum is therefore the initial concentration, not a delayed peak. For 0<<1, the algebraic tail makes the infinite-horizon area under the concentration curve diverge. In the illustrative coupled simulations, the apparent late tumour reduction at small is consequently driven by persistent modelled drug exposure and must not be read as evidence of superior therapy. All parameters are dimensionless and illustrative; no patient data are fitted. Source code and generated tables are included to make every numerical statement reproducible.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

A major transition in the evolution of eukaryotic cells was phagocytosis: the ability of one cell to engulf another. In some cases, the engulfed cell persists, forming an endosymbiosis which is often mutually beneficial. However, it is unclear whether the interaction between the ancestors was also cooperative. There are clear benefits to developing the ability to eat a competitor, but can it ever be beneficial to eat a cooperator? We use a dynamical systems framework to identify the ecological conditions under which phagocytosis can evolve from a system of cooperating microbes.Beginning with a system of two cooperating species, we use invasion analysis to identify the parameter regimes under which a phagocytosing mutant can invade and the subsequent population dynamics. We identify a critical threshold that determines whether an invading predator will coexist with the cooperative community or replace members. Repeatedly invading with different predators, we find that the system saturates when there are three coexisting species, though this community does not have a single optimal composition. However, if we incorporate additional costs, we find there is a single evolutionarily stable community of two cooperators and a pure predator. Ultimately, we show that it is possible for phagocytosis to evolve from a cooperative community and, moreover, coexist with the original cooperating species.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Buruli ulcer (BU) is a neglected tropical skin disease caused by infection with Mycobacterium ulcerans. The mode of its transmission into the skin has not been fully clarified, but it is believed to involve mechanical transfer, perhaps by insect bites. M. ulcerans secretes a toxin called mycolactone which causes necrosis of the skin and formation of large skin lesions. Understanding the development of these skin lesions over time is crucial to understand the pathogenesis of BU. We formulate a mathematical model to describe the space-time dynamics of M. ulcerans bacteria, mycolactone and skin cell densities as well as fibrin deposition. We investigate the relative contributions to ulcer formation of direct cytotoxicity of mycolactone and ischaemia following the loss of coagulation control. Our results are consistent with existing biological evidence that higher mycolactone concentrations accelerate cell death. Furthermore, the model indicates that the rate of wound expansion is primarily controlled by the diffusion and net birth rate of bacteria, which sets the overall scaling for the speed of wound expansion, along with two dimensionless constants that define the shape of the trajectory of wound edge position against time, namely the ratio of the mycolactone diffusion coefficient to the bacterial diffusion coefficient, and the ratio of the mycolactone degradation rate to the bacterial net birth rate. When mycolactone diffuses significantly faster than bacteria spread or when mycolactone degrades more slowly than bacteria reproduce, the initial development of the wound is faster. In contrast, the long-term wound expansion speed appears to be relatively unchanged across a wide range of the dimensionless parameters.

Species at the highest trophic levels are generally the least resilient to perturbations, yet the factors that determine their persistence remain poorly understood. We ask how the intrinsic self-regulation of a basal host governs the persistence of the trophic levels above it. To address this question, we use a minimal discrete-time model of a host-parasitoid-hyperparasitoid system motivated by the biology of trigonalid wasps, obligate hyperparasitoids that complete development only when their caterpillar host is also attacked by a primary parasitoid. We compare two forms of host self-regulation in the absence of natural enemies: scramble (Ricker) fecundity, appropriate when the parasitoid is the host's principal enemy and which recovers the classical model of Beddington and Hammond as a special case; and contest (Beverton-Holt) fecundity, appropriate when the host is held in check by many additional enemies. Using numerical continuation, Lyapunov exponents, and direct computation of the model's trajectories, we find that the form of host regulation is decisive for persistence at the top trophic level. Under scramble host dynamics, coexistence is extremely fragile with respect to both parameter values and initial conditions. Under contest host dynamics, coexistence is robust across broad parameter ranges and loses stability only gradually, through a supercritical Neimark--Sacker bifurcation and a subsequent torus-doubling route to chaos. Therefore, in at least some trophic networks, the resilience of the highest levels can be determined at the base of the network.

Bacterial mutation rates are traditionally inferred from phenotypic data using fluctuation assays. Mutation rate computation from these assays relies on a mathematical model describing the emergence of mutants during population growth. However, this standard model relies on restrictive assumptions that are often transgressed in real biological setups. Key assumptions include equal growth rate of the wild-type and the mutant population, no cell death during the growth phase, and no post-growth sampling. Although several refined mathematical models have been proposed over the last decades to circumvent some of these assumptions, none can fully account for the complex ecological scenarios frequently observed in the lab or in the wild, such as the combination of cell death and an unknown fitness effect of the mutation, or non-constant death rates. In this work, we propose to replace these standard mathematical models by a stochastic simulation framework, which can circumvent all these assumptions. We provide a fast and accurate implementation of this stochastic model, and combine it with an Approximate Bayesian Computation scheme to perform simulation-based, likelihood-free parameter inference. Our approach can simultaneously infer mutation rate and fitness effect of the mutation while accounting for death and partial plating, and could be extended to arbitrary complex conditions that can be efficiently simulated.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Historical plague mortality records provide unusually detailed temporal records of epidemic dynamics. However, the ecological processes underlying those dynamics are rarely observed. This creates a fundamental challenge for inference into ecological mechanisms, as multiple combinations of unobserved ecological parameters and initial conditions may produce similar mortality trajectories. We investigated this problem by examining the geometry of parameter compensation for models of the Barcelona plague outbreak of 1490. We fit a mechanistic rat-flea-human transmission model to daily human mortality data under three values of flea search efficiency, then characterized the resulting joint posterior distributions. We then used a phenomenological model in which human infection was driven by a time-varying external force of infection, allowing us to separate parameter dependencies arising from deterministic structure. The mechanistic model reproduced the broad temporal structure of the observed epidemic, but exhibited strong linear and nonlinear dependencies among transmission parameters and initial rat and flea conditions. Similar dependencies persisted in the phenomenological model, despite its substantially reduced complexity. However, fixing the timing of the external hazard substantially reduced these dependencies. These results demonstrate that detailed mortality records can constrain the temporal pattern of an epidemic without uniquely identifying the mechanisms that generated it. Rather than simply adding additional mortality records, improving mechanistic inference from historical epidemics may require recovering independent observations of upstream ecological processes, such as timing of rodent or flea dynamics as well as environmental drivers. Mechanistic models can therefore serve not only to reconstruct historical epidemics but also to identify which missing observations would most improve inference.

In this study, we develop a mathematical model of the hypothalamic pituitary adrenal (HPA) axis and systematically investigate the model complexity required to reproduce observed hormone dynamics. The model is fitted to ACTH and cortisol measurements from five healthy volunteers reported in the literature. We show that simultaneous reproduction of both hormone profiles requires a delay differential equation model with at least five state variables among the model architectures considered. Following fitting the volunteer data, we investigate the response of the model to sustained stress using bifurcation analysis. The coupled system exhibits bifurcations over a range of stress intensities. This analysis reveals bifurcations whose critical stress thresholds depend strongly on the modulatory state of the system. These results provide a mathematical framework for investigating transitions between physiological and dysregulated HPA axis dynamics, which may be related to states of health and disease.

Oncogene amplification on circular extrachromosomal DNA (ecDNA) has been linked to poor prognosis and higher treatment resistance in multiple types of human cancer. ecDNA are mobile genetic elements lacking centromeres that are partitioned unevenly into daughter cells at mitosis. While random segregation of ecDNA contributes to gene-copy-number heterogeneity among tumour cells, how ecDNA contribute to phenotypic heterogeneity, via mRNA and proteins, and how this affects targeted therapy outcomes is still not understood. In fact, cancer cell populations with ecDNA show remarkable heterogeneity in mRNA concentrations, which is not explained by gene-copy-number heterogeneity alone. In this paper, we use a stochastic model of gene transcription to mathematically model mRNA production and decay in cancer cells with ecDNA, extending a method that has been used to study phenotypic heterogeneity in eurokaryotic and prokaryotic cells. Using our model, we explore the space of possible mRNA distributions from a given ecDNA distribution, under different model parameter regimes. Our model is able to produce realistic mRNA distributions and we use it to estimate transcriptional parameters from single-cell measurements of ecDNA and mRNA counts in cell-line data. We discover a remarkable similarity in estimated parameters across cell lines and amplification types, indicating similar processes underlying transcription in ecDNA and chromosomal amplifications. We also model genotoxic therapy and use our model to predict the mRNA composition of persisting cells, showing the usefulness of such a modelling framework to understanding the origins of treatment resistance in ecDNA-containing cancers.

The mechanisms by which organisms extract relevant information from complex sensory signals and use it for decision-making constitute a fundamental issue in biological and cognitive sciences. In social animal societies, individual decision-making is profoundly shaped by information provided by conspecifics. Collective patch-foraging experiments in the laboratory provide a controlled setting in which social information use can be quantified. Here, we consider a go/no-go task in which groups choose between two patches differing in food reward probability. We use agent-based simulations, underpinned by an augmented collective drift-diffusion model, to investigate alternative mechanisms for representing and integrating social information. We consider two representations, continuous (counting representation) and discrete (pulsatile representation), and two integration mechanisms, biasing the decision threshold (threshold modulation) and modifying the accumulated belief (belief modulation), yielding four distinct social models, in addition to a non-interacting model. We first characterize the collective dynamics generated by these models across cognitive parameters and experimental conditions. The temporal dynamics of group accuracy provide informative signatures of the underlying mechanisms. We then assess model selection and parameter identifiability using Bayesian inference and Wasserstein distance minimization, applied to group-accuracy and departure-time distributions. Bayesian inference using distributions of group accuracy provides the most reliable identification across the conditions considered. Importantly, model selection and parameter identifiability are associated with different aspects of collective behavior: model selection is closely linked to the presence of oscillations in the temporal dynamics, whereas parameter identifiability is more closely related to the global accuracy value. Thus, the information available for distinguishing the underlying cognitive mechanisms is not necessarily the same as that required to recover their parameters. Our results further show that identification depends not only on model structure but also on the experimental conditions. Overall, this study provides practical guidelines for identifying socio-cognitive mechanisms in collective foraging experiments.

Cancer progression often requires the accumulation of multiple driver mutations, yet the same set of drivers may be acquired in different orders. How these alternative mutation-order pathways jointly shape tumor clonal structure remains unclear. We develop a multitype branching-process model in which malignant transformation requires two driver mutations and distinguish malignant cells both by their mutation order and by the independent transformation event that founded their clone. Under a successive exponential approximation, we establish point-process limits for the pathway-specific clone-size processes and derive a closed-form expression for the limiting expected Simpson's index of the combined malignant population. When the two mutation orders produce malignant cells with the same net growth rate, the limiting index separates into effective pathway weights, determined by mutation rates and birth-death dynamics at preceding stages, and within-pathway concentration terms, determined by intermediate-to-malignant growth-rate ratios. This decomposition shows that a driver's effect on heterogeneity depends critically on when it is acquired. A strong driver acquired early expands the intermediate lineage and increases the supply of independent malignant founders, whereas the same driver acquired at the final transition strengthens the growth and age advantage of early-founded malignant clones. Under additive fitness effects, these counteracting mechanisms can produce a non-monotone relationship between selective advantage and clonal concentration. We further show that threshold-like non-additive fitness effects can generate highly concentrated malignant populations, while order-dependent terminal fitness leads the faster-growing pathway to dominate asymptotically. Together, these results reveal how mutation order, mutational accessibility, selection, and epistasis jointly determine lineage-level intratumor heterogeneity.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Many pathogenic and microbial species undergo complex life cycles wherein they experience drastic bottlenecks during invasion into a host and rapid growth within hosts. The effects of such a complex population history on patterns of genetic variation is not currently understood. We model a pathogen population as a metapopulation where each deme represents a host and incorporate transmission dynamics as well as the effects of transmission bottlenecks. We employ a coalescent framework, obtaining analytical expressions for pairwise times to coalescence and genetic differentiation within and between hosts. We find that recurrent bottlenecks rescale the coalescent process within hosts to the size of the bottleneck, reducing within-host variation. In addition, the pairwise distribution of time to coalescence under our model deviates substantially from the Wright-Fisher process when the number of hosts (d) is of similar order to the size of the bottleneck (k), i.e., k [~] d and when m <<1/(2k) +1/(2d), where m is the migration rate. Otherwise, the pairwise process under our model converges in distribution to the Wright-Fisher process. Coalescent simulations show that under the same conditions the site frequency spectrum (SFS) under our model deviates from that under the Kingman coalescent and is drastically skewed towards singletons. Our results suggest that the inference of selection or demographic history from pathogen population genetics data may be systematically biased by the mechanics of pathogen reproduction, indicating a need for further theoretical work that can inform inference.

Biologists often wish to explain why one trait is more common than its alternatives among species in a clade: is it due to asymmetries in diversification or character evolution processes? State-dependent speciation-extinction (SSE) models provide a flexible framework for jointly modeling such processes. Here, we derive mathematical results that define functional relationships (``rulesets") among rate parameters, establish a mapping between rulesets and evolutionary scenarios, and characterize how these scenarios give rise to different long-term tip state distributions. In particular, we show that multiple evolutionary scenarios can produce identical stationary distributions of tip state frequencies. Performing inference on simulated data from extant-only trees, we find that not all scenarios are equally favored for explaining datasets constructed solely from extant species. This bias could arise from the interplay between mixing rates, elapsed time, and the conditioning inherent in tree likelihoods. As a proof of concept, we apply our framework to the binary SSE (BiSSE) model across multiple Squamate clades and conduct a literature survey of previous BiSSE analyses on Angiosperm phylogeny. We show that, to some extent, the preference for certain scenarios observed in simulations is also reflected in analyses of empirical trees. This study provides a theoretical ground to directly link data-generating process and observed tip state frequencies under an SSE framework.

The Yampolsky-Stoltzfus (YS) model provided the first formal demonstration that a bias in the introduction of new mutations can shape the direction of adaptive evolution, even when mutation rates are small relative to selection coefficients. Despite its foundational role, the model has been characterized almost entirely by simulation and closed-form solutions for its central quantity, the probability PA that a mutationally favored but less-fit genotype fixes ahead of a rarer, fitter alternative, have been available only in the limiting regimes of vanishing mutation supply and infinite population size. Here we derive closed-form approximations for PA across the intermediate regime of finite population size and appreciable mutation supply that lies between these two limits. Under strong selection between the competing mutants, we combine an establishment-probability argument with clonal interference theory to obtain an explicit expression for PA that accounts for both the ongoing input of the mutationally favored mutation and the frequency-dependent establishment probability of its competitor. Under weak selection between competing mutants, we derive a complementary expression based on the fixation of a mixed population under Kimura's classical formula. Both expressions agree closely with direct simulation of the underlying Wright-Fisher process, recover the known limiting cases, and remain consistent with one another over a substantial intermediate range of parameters. These results are expected to provide a basis for assessing how the strength of mutation bias affects adaptation in empirical data.

Game theory in biology started as an attempt to model animal contests by making assumptions about the fitness costs individuals pay in aggressive interactions. Among the different approaches, a notable one is to construct mechanistic models, using assumptions about behaviours and cognitive processes. Reinforcement learning is an important cognitive process, and we use it here to model pairwise and coalitional contests. We study situations where territorial neighbours become so-called dear enemies. A defending individual can get help from a neighbour through a defender-neighbour coalition against a challenger attempting territory takeover. In our model, coalition members have an advantage in contests in terms of costs of aggression, with the challenger being exposed to aggression from both members and each member only receiving part of the aggression from the challenger. We find that notable costs of border conflicts between neighbours, together with substantial advantages for coalitions, favour intervention. Neighbours intervene when they can help weaker defenders, in this way avoiding costs of border renegotiation with a new and potentially stronger territory owner. We also introduce different forms of perceived costs into the model, which we show correspond to the range from pure mutual assessment to partial self-assessment.